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THEORETICAL  ANALYSIS  OF  MICROWAVE  AND  MILLIMETER  WAVE  INTEGRATED 
CIRCUITS  BASED  ON  MAGNETIC  FILMS 

Under  the  sponsorship  of  the  ONR  Contract  Contract  N00014-89-J-1019  we  have 
published  29  refereed  journal  and  conference  papers. 

A  macroscopic  model  is  proposed  to  explain  nonlinear  electromagnetic  phenomena  in 
superconductors.  Nonlinear  constitutive  relations  for  electromagnetic  problems  are  derived 
by  modifying  the  linear  London’s  equations.  The  superelectron  number  density  na  is  a 
function  of  the  applied  current  density.  The  critical  current  density  Jc  is  derived  classically 
from  a  critical  energy  Ec.  For  temperature  T  /  0,  the  concept  of  critical  current  Jc 
does  not  imply  an  abrupt  transition  of  the  whole  sample  from  a  superconducting  state 
to  a  normal  state  when  J  >  Jc.  A  rather  smooth  variation  of  na(J)  is  shown  instead. 
The  relation,  na(J),  is  derived  from  the  Maxwellian  distribution  of  electron  velocities 
at  a  certain  temperature  T  and  a  certain  macroscopic  current  density  J.  Agreement 
has  also  been  found  between  this  na(J,T)  model  and  the  temperature  dependance  of 
na  in  the  two-fluid  model.  The  nonlinear  conductivities  <r3( J)  and  <rn(J )  are  obtained 

from  the  London’s  equation  and  the  na(J)  function.  Nonlinear  resistance  i2(/),  kinetic 

!  ""  ’ 

inductance  Ljg(I).  and  surface  impedance  Za(I)  in  thin  wire,  slab,  and  strip  geometries 
of  superconductors  are  calculated.  A  general  scheme  of  solving  nonlinear  electromagnetic 
problems  in  superconductors  is  proposed.  A  good  agreement  between  the  theory  and 
experiments  has  been  found. 

The  Riccati  differential  equation  for  reflection  coefficients  in  one-dimensional  in¬ 
homogeneous  media  is  applied  to  the  electromagnetic  inverse  scattering  problem.  Two 

h~. 

types  of  Riccati  equation  in  literature,  Schelkunoff ’s  and  Redheffer’s,  are  derived  and  dis¬ 
tinguished.  Based  on  inverting  Redheffer’s  Riccati  equations,  both  linear  and  non-linear 
inversion  formulae  are  proposed.  These  renormalized  perturbation  formulae  reconstruct 
the  dielectric  profile  from  the  reflection  coefficient  at  the  surface  of  the  medium.  Exist¬ 
ing  inversion  formulae  (including  the  Born  approximation)  which  were  obtained  from  the 
Green’s  function  approach  and  the  Gel’fand-Levitan-Marchenko  (GLM)  theory  are  now 
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derived  from  the  Riccati  equation.  Four  inversion  schemes  based  on  inverting  the  Riccati 
equation,  linearized  Rcdheffer’s  Riccati  equation,  linearized  Schelkunoff’s  Riccati  equation, 
non-linear  Redheffer’s  Riccati  equation  and  non-linear  Schelkunoff’s  Riccati  equation  ap¬ 
proaches,  are  used  to  invert  several  dieletric  profiles.  Comparison  and  summary  of  these 
methods  are  given.  All  these  methods  give  higher  order  inversion  results  than  the  first-order 
Born  approximation.  The  inversion  is  performed  on  band-limited  reflection  coefficients  in 
frequency  domain.  An  inverse  Liouville  transform  is  introduced  to  rigorously  recover  the 
geometric  lengths  from  the  stretched  coordinates  in  the  inversion  procedure. 

In  general,  the  inverse  methods  can  be  summarized  in  three  categories.  The  first 
one  is  a  data  base  approach.  Whenever  a  measurement  result  comes  in,  it  will  be  compared 
with  all  the  pre-stored  patterns  in  a  library  to  identify  the  possible  targets.  In  industry,  it  is 
called  data  interpretation  which  strongly  relies  on  previous  experience.  The  accumulation 
of  data  in  the  library  is  a  learning  process.  The  advantage  of  this  method  is  that  the 
unknowns  can  be  identified  exactly  if  such  target  has  been  seen  before  so  that  the  cause 
and  the  result  are  directly  connected.  The  disadvantage  is  that  if  the  data  base  is  not  big 
enough,  the  target  can  not  be  identified,  or  if  the  data  base  is  too  big,  it  will  take  too 
much  time  to  search  for  and  match  the  right  pattern.  The  second  method  is  the  iterative 
approach.  A  forward  model  of  deducing  results  from  cause  has  been  established  before  the 
inversion.  When  a  measurement  is  obtained,  a  set  of  guessed  values  for  the  unknowns  will 
be  put  in  the  forward  model  as  the  possible  causes.  If  the  predicted  result  of  the  forward 
model  does  not  match  the  measurement,  an  adjusted  set  of  parameters  will  be  used  as 
a  new  guess.  The  whole  procedure  will  be  repeated  until  the  match  is  found.  During 
iterations,  how  to  correct  the  error  to  reach  fastest  the  convergent  value  is  very  important. 
This  procedure  is  also  called  optimization  or  error  minimization,  which  by  itself  is  an  active 
research  area.  Iterative  methods  are  particularly  useful  when  the  solution  of  an  integral 
equation  can  not  be  found  in  an  explicit  form.  As  an  example  of  the  iterative  methods,  the 
Born  approximation  approximates  the  forward  model  in  a  linearized  form,  then  iterates 
towards  convergence.  The  advantage  of  the  iterative  approach  is  that  an  accurate  final 
solution  can  be  found  at  a  relative  fast  (compare  to  data  base  search)  speed.  The  method 
is  also  flexible  for  different  kinds  of  targets,  unlike  the  data  base  search  where  the  target 
has  to  have  been  seen  before.  The  disadvantage  is  that  the  convergence  of  iteration  is 
not  always  guaranteed.  There  are  usually  more  than  one  minimum  of  the  cost  function. 


3 


The  initial  guess  is  sometimes  so  critical  to  assure  the  right  convergent  result  that  a  Priori 
information  needs  to  be  used.  Non-uniqueness  is  a  serious  problem  in  all  inversion  methods. 
Sufficient  number  of  measurements  are  helpful  of  reducing  the  degree  of  non-uniqueness. 
Formulation  of  the  forward  model  for  the  final  inversion  equation  determines  how  efficient 
the  inversion  scheme  is,  which  sometimes  also  affect  the  degree  of  non-uniqueness. 

A  recently  developed  inversion  method  referred  to  as  the  renormalized  Source- Type 
Integral  Equation  (STIE)  approach  solves  the  integral  equation  derived  from  the  Green’s 
function  without  the  linear  approximation.  The  STIE  approach  formulates  an  exact  for¬ 
ward  model,  hence  is  not  restricted  to  low  contrast  profiles  or  weak  scattering.  The  STIE 
method  has  been  applied  to  profile  inversion  problems  of  the  soil  moisture  and  oil  for¬ 
mation  in  boreholes,  where  the  medium  properties  are  described  by  the  permittivity  and 
conductivity  distributions.  The  unknown  profiles  are  inverted  from  the  electromagnetic 
measurements  at  remote  observation  points.  The  third  category  of  inverse  methods  is  the 
explicit  solution  approach.  When  one  solves  the  inversion  equation,  closed-form  solutions 
may  be  obtained  if  the  problem  is  formulated  in  certain  ways  and  the  data  is  in  certain 
forms.  Examples  of  this  method  include  the  Gel’fand-Levitan-Marchenko  (GLM)  theory 
when  rational  reflection  coefficient  is  used  and  methods  of  inverting  the  Riccati  equation. 
The  advantage  is  obviously  that  the  convergence  is  guaranteed  and  the  speed  is  excellent. 
Unfortunately,  only  a  very  few  practical  problems  can  be  solved  by  using  this  method.  It 
is  usually  used  as  a  check  for  other  inversion  methods  or  for  mathematical  studies.  The 
STIE  approach  is  extended  to  the  case  of  a  general  background  medium.  The  formulation 
of  the  inversion  equation  is  exact.  Numerical  optimization  methods  are  employed  in  the 
solution  of  the  inversion  equation.  The  inversion  equation  also  has  an  explicit  dependence 
on  the  unknowns  to  be  inverted  for;  this  allows  one  to  compute  the  derivative  (including 
higher  orders)  of  the  response  with  respect  to  these  unknowns  in  a  closed  form.  Therefore, 
higher-order  convergence  can  be  achieved  without  increasing  much  of  the  computational 
time.  By  pre-storing  the  elements  of  the  inversion  equation  (which  depend  only  on  the 
background  medium  and  are  independent  of  the  unknown  profile),  the  method  does  not 
require  the  solution  to  the  full  forward  problem  repeatedly  as  in  the  case  of  the  Distorted 
Born  approach  and  is  therefore  faster  in  implementation.  An  algorithm  is  formulated 
for  simultaneously  inverting  the  unknown  permittivity  and  conductivity  profiles.  As  an 
example,  a  soil  moisture  profile  is  inverted  from  measurements  above  the  ground. 


A  full  modal  analysis  is  used  to  study  the  dispersion  characteristics  of  microstrip 
lines  periodically  loaded  with  crossing  strips  in  a  stratified  uni  axially  anisotropic  medium. 
Dyadic  Green’s  functions  in  the  spectral  domain  for  the  multilayered  medium  in  conjunc¬ 
tion  with  the  vector  Fourier  transform  (VFT)  are  used  to  formulate  a  coupled  set  of  vector 
integral  equations  for  the  current  distribution  on  the  signal  line  and  the  crossing  strips. 
Galerkin’s  procedure  is  applied  to  derive  the  eigenvalue  equation  for  the  propagation  con¬ 
stant.  The  effect  of  anisotropy  for  both  open  and  shielded  stuctures  on  the  stopband 
properties  is  investigated. 

The  input  impedance  of  a  microstrip  antenna  consisting  of  two  circular  microstrip 
disks  in  a  stacked  configuration  driven  by  a  coaxial  probe  is  investigated.  A  rigorous  anal¬ 
ysis  is  performed  using  a  dyadic  Green’s  function  formulation  where  the  mixed  boundary 
value  problem  is  reduced  to  a  set  of  coupled  vector  integral  equations  using  the  vector 
Hankel  transform.  Galerkin’s  method  is  employed  in  the  spectral  domain  where  two  sets 
of  disk  current  expansions  are  used.  One  set  is  based  on  the  complete  set  of  orthogonal 
modes  of  the  magnetic  cavity,  and  the  other  employs  Chebyshev  polynomials  with  the 
proper  edge  condition  for  the  disk  currents.  An  additional  term  is  added  to  the  disk  cur¬ 
rent  expansion  to  properly  model  the  current  in  the  vicinity  of  the  probe/disk  junction. 
The  input  impedance  of  the  stacked  microstrip  antenna  including  the  probe  self-impedance 
is  calculated  as  a  function  of  the  layered  parameters  and  the  ratio  of  the  two  disk  radii. 
Disk  current  distributions  and  radiation  patterns  are  also  presented.  The  calculate  results 
are  compared  with  experimental  data  and  shown  to  be  in  good  agreement. 

The  coupled-wave  theory  is  generalized  to  analyze  the  diffraction  of  waves  by  chiral 
gratings  for  arbitrary  angles  of  incidence  and  polarizations.  Numerical  results  for  the 
Stokes  parameters  of  diffracted  Floquet  modes  versus  the  thickness  of  chiral  gratings  with 
various  chiralities  are  calculated.  Both  horizontal  and  vertical  incidences  sire  considered  for 
illustration.  The  diffracted  waves  from  chiral  gratings  are  in  general  elliptically  polarized; 
and  in  some  particular  instances,  it  is  possible  for  chiral  gratings  to  convert  a  linearly 
polarized  incident  field  into  two  nearly  circularly  polarized  Floquet  modes  propagating  in 
different  directions. 
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A  general  spectral  domain  formulation  to  the  problem  of  radiation  of  arbitrary 
distribution  of  sources  embedded  in  a  horizontally  stratified  arbitrary  magnetized  linear 
plasma  is  presented.  The  fields  are  obtained  in  terms  of  electric  and  magnetic  type  dyadic 
Green’s  functions.  The  formulation  is  considerably  simplified  by  using  the  kDB  system 
of  coordinates  in  conjunction  with  the  Fourier  transform.  The  distributional  singular 
behavior  of  the  various  dyadic  Green’s  functions  in  the  source  region  is  investigated  and 
taken  into  account  by  extracting  the  delta  function  singularities.  Finally,  the  fields  in  any 
arbitrary  layer  are  obtained  in  terms  of  appropriately  defined  global  upward  and  downward 
reflection  and  transmission  matrices. 

We  have  investigated  a  method  for  the  calculation  of  the  current  distribution,  re¬ 
sistance,  and  inductance  matrices  for  a  system  of  coupled  superconducting  transmission 
lines  having  finite  rectangular  cross  section.  These  calculation  allow  accurate  characteri¬ 
zation  of  both  high-Tc  and  1ow-jTc  superconducting  strip  transmission  lines.  For  a  single 
stripline  geometry  with  finite  ground  planes,  the  current  distribution,  resistance,  induc¬ 
tance,  and  kinetic  inductance  are  calculated  as  a  function  of  the  penetration  depth  for 
various  film  thickness.  These  calculations  are  then  used  to  determine  the  penetration 
depth  for  Nb,NbN,  and  Y  BaiCu^O’i-x  superconducting  thin  films  from  the  measured 
temperature  dependence  of  the  resonant  frequency  of  a  stripline  resonator.  The  calcula¬ 
tions  are  also  used  to  convert  measured  temperature  dependence  of  the  quality  factor  to 
the  intrinsic  surface  resistance  as  a  function  of  temperature  for  a  Nb  stripline  resonator. 

Proximity-print  x-ray  lithography  is  commonly  performed  with  gold  or  tungsten 
structures  of  sizes  down  to  30  nm  wide  and  50-800  nm  tall  which  are  patterned  onto  the 
surface  of  a  thin,  x-ray  transparent  membrane.  X-rays  in  the  wavelength  range  of  0.5-5  nm 
sure  used  for  replication  with  mask-substrate  gaps  ranging  from  zero  (contact  print)  up  to 
20  /xm  or  more.  The  resolution  o  f  this  method  (minimum  achievable  linewidth)  is  limited 
predominantly  by  the  diffraction  of  the  x-rays  around  these  structures  and  the  spreading 
of  the  diffracted  waves  into  the  0-20  fim  gap.  Work  to  date  has  assumed  that  scalar 
diffraction  theory  is  applicable — as  calculated,  for  example,  by  the  Rayleigh-Sommerfeld 
formulation — and  that  Kirchhoff  boundary  conditions  can  be  applied.  Kirchhoff  boundary 
conditions  assume  that  the  fields  are  constant  in  the  region  between  the  absorbers,  and 
also  (a  different)  constant  in  the  region  just  under  the  absorbers,  and  that  there  are  no 
fringing  fields.  In  this  report  we  explore  the  validity  of  this  assumption  for  the  case 
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of  30  nm-wide  by  30-100  nm-tall  gold  absorbers  with  4.5  nm  (Cg)  x-rays.  Because  of 
computational  time  limitations,  the  shorter  wavelength  and  larger  absorber  cases  are  not 
currently  possible.)  Because  the  absorber  is  only  7  wavelengths  wide  and  7-20  wavelengths 
high,  strong  diffractive  effects  are  expected.  The  finite- difference  time-domain  (FD-TD) 
technique  was  used  on  a  Cray- 2  supercomputer  to  predict  the  fields  diffracted  by  the  gold 
absorbers.  In  applying  the  FD-TD  technique,  Maxwell’s  equations  are  discretized  in  space 
and  time  on  a  uniform  rectangular  grid.  A  second-order  absorbing  boundary  condition  is 
applied  at  the  outer  boundary  of  the  computational  domain  in  order  to  simulate  unbounded 
space.  The  results  indicate  that  strong  fringing  fields  exist  in  the  shadow  region  of  the 
absorber,  and  hence  Kirchhoff  boundary  conditions  are  not  accurate  in  this  regime. 

Because  the  effects  of  diffraction  during  proximity-print  x-ray  lithography  are  of 
critical  importance,  a  number  of  previous  researchers  have  attempted  to  calculate  the 
diffraction  patterns  and  minimum  achievable  feature  sizes  as  a  function  of  wavelength  and 
gap.  Work  to  date  has  assumed  that  scalar  diffraction  theory  is  applicable-as  calculated,  for 
example,  by  the  Rayleigh-Sommerfeldformulation-andthat  Kirchhoff  boundary  conditions 
can  be  applied.  Kirchhoff  boundary  conditions  assume  that  the  fields  (amplitude  and 
phase)  are  constant  in  the  open  regions  between  absorbers,  and  a  different  constant  in 
regions  just  under  the  absorbers  (i.e.,  that  there  are  no  fringing  fields).  An  x-ray  absorber 
is,  however,  best  described  as  a  lossy  dielectric  that  is  tens  or  hundreds  of  wavelengths 
tall,  and  hence  Kirchhoff  boundary  conditions  are  unsuitable.  In  this  report  we  use  two 
numerical  techniques  to  calculate  (on  a  Cray  2  supercomputer)  accurate  diffracted  fields 
from  gold  absorbers  for  two  cases:  a  30  nm-wide  line  at  X  =  4.5  nm,  and  a  100  nm-wide 
line  at  A  =  1.3  nm.  We  show  that  the  use  of  Kirchhoff  boundary  conditions  introduces 
unphysically  high  spatial  frequencies  into  the  diffracted  fields.  The  suppression  of  these 
frequencies-which  occurs  naturally  without  the  need  to  introduce  an  extended  source  or 
broad  spectrum-improves  exposure  latitude  for  mask  features  near  0.1  /im  and  below. 
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In  order  to  understand  the  physical  meaning  of  rational  reflection  coefficients  in 
one-dimensional  inverse  scattering  theory  for  optical  waveguide  design,  we  have  studied 
the  relation  between  the  poles  of  the  transverse  reflection  coefficient  and  the  modes  in 
inhomogeneous  dielectrics.  By  using  a  stratified  medium  model  it  is  shown  that  these 
poles  of  the  reflection  coefficient  have  a  one-to-one  correspondence  to  the  discrete  modes, 
which  are  the  guided  and  leaky  modes.  The  radiation  modes  have  continuous  real  values  of 
transverse  wave  numbers  and  lire  not  represented  by  the  poles  of  the  reflection  coefficient. 
Based  on  these  results,  applications  of  the  Gel’fand-Levitan-Marchenko  theory  to  optical 
waveguide  synthesis  with  the  rational  function  representation  of  the  transverse  reflection 
coefficient  are  discussed. 

We  developed  an  inversion  algorithm  based  on  a  recently  developed  inversion  method 
referred  to  as  the  Renormalized  Source-Type  Integral  Equation  approach.  The  objective 
of  this  method  is  to  overcome  some  of  the  limitations  and  difficulties  of  the  iterative  Born 
technique.  It  recasts  the  inversion,  which  is  nonlinear  in  nature,  in  terms  of  the  solution 
of  a  set  of  linear  equations;  however,  the  final  inversion  equation  is  still  nonlinear.  The 
derived  inversion  equation  is  an  exact  equation  which  sums  up  the  iterative  Neuman  (or 
Born)  series  in  a  closed  form  and;  thus,  is  a  valid  representation  even  in  the  case  when 
the  Born  series  diverges;  hence,  the  name  Renormalized  Source-Type  Integral  Equation 
Approach. 

The  scattering  and  receiving  characteristics  of  a  probe-fed  stacked  circular  mi¬ 
crostrip  antenna,  both  as  an  isolated  element  and  in  an  infinite  array,  are  investigated. 

The  receiving  case,  where  the  antenna  is  loaded  with  impedance  Z^,  is  solved  by 
superposition,  decomposing  the  problem  into  the  scattering  case  with  Z ^  =  0  and  the 
transmitting  case.  In  the  scattering  case,  the  coaxial  probe  is  short-circuited  to  the  ground 
plane  and  the  induced  probe  current  I\  due  to  an  incident  plane  wave  excitation  is  de¬ 
termined.  In  the  transmitting  case,  a  voltage  V  is  applied  to  the  base  of  the  probe  and 
the  input  impedance  Z'm  is  solved  for,  giving  a  relationship  between  the  applied  voltage 
V  and  the  transmitting  probe  current  I2-  With  the  knowledge  of  I\  and  ^in>  for  a  given 
load  impedance  the  total  probe  current,  I  =  I\  -f  I2,  and  the  received  power  are 
determined. 
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The  scattering  and  transmitting  problems  are  solved  rigorously  using  a  dyadic 
Green’s  function  formulation  where  the  mixed  boundary  value  problem  is  reduced  to  a  set 
of  coupled  vector  integral  equations  for  the  unknown  disk  and  probe  currents.  Galer kin’s 
method  is  employed  in  the  spectral  domain  where  the  disk  current  distributions  are  ex¬ 
panded  in  terms  of  the  complete  set  of  transverse  magnetic  (TM)  and  transverse  electric 
(TE)  modes  of  a  cylindrical  resonant  cavity  with  magnetic  side  walls.  An  additional  term 
is  added  to  the  disk  current  expansion  to  properly  model  the  singular  behavior  of  the 
current  in  the  vicinity  of  the  probe,  to  ensure  continuity  of  the  current  at  the  probe/disk 
junction,  and  to  speed  up  the  convergence  of  the  solution. 

The  radar  cross  section  (RCS)  of  a  single  stacked  microstrip  antenna  is  calculated 
for  both  the  open  and  short-circuited  cases.  For  an  infinite  array  of  phased  elements,  the 
reflection  coefficient  seen  at  the  input  of  the  antenna  and  the  received  power  are  calculated. 

The  complex  resonant  frequencies  of  the  open  structure  of  a  microstrip  antenna 
consisting  of  two  circular  microstrip  disks  in  a  three  layer  stacked  configuration  have  been 
rigorously  calculated  as  a  function  of  the  layered  parameters  and  the  ratio  of  the  radii 
of  the  two  disks.  Using  a  dyadic  Green’s  function  formulation  for  horizontally  stratified 
media  and  the  vector  Hankel  transform,  the  mixed  boundary  value  problem  is  reduced  to 
a  set  of  coupled  vector  integral  equations.  Employing  Galerkin’s  method  in  the  spectral 
domain,  the  complex  resonant  frequencies  are  calculated  and  convergence  of  the  results  is 
demonstrated.  It  is  shown  that  for  each  mode,  the  stacked  circular  microstrip  structure  has 
dual  resonant  frequencies  which  are  associated  with  the  two  coupled  constitutive  resonators 
of  the  structure  and  which  are  a  function  of  the  mutual  coupling  between  them.  This  mutal 
coupling  depends  on  the  geometrical  configuration  of  the  stacked  structure,  the  layered 
parameters,  and  the  disk  radii.  The  maximum  coupling  effect  occurs  where  the  real  parts 
of  the  resonant  frequencies  of  the  constitutive  resonators  are  approximately  equal,  where 
the  behavior  of  the  resonances  in  this  region  is  a  function  of  the  coupling.  The  dual 
frequency  behavior  of  the  stacked  microstrip  structure,  easily  controlled  by  varying  the 
parameters  of  layer  2  and  disk  radii  ratio,  given  fixed  parameters  for  layer  1  and  layer  3, 
may  be  used  to  broaden  the  bandwidth  or  provide  for  dual  frequency  use  of  the  antenna. 
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We  rigorously  analyze  the  radiation  problem  of  a  circular  patch  which  is  center  fed 
by  a  coaxial-line  driven  probe  over  a  ground  plane  and  situated  in  an  arbitrary  layered 
medium.  The  current  distribution  on  both  the  patch  and  the  probe  is  rigorously  formu¬ 
lated  using  a  planar  stratified  medium  approach.  A  set  of  three  coupled  integral  equation 
is  derived  which  governs  the  axial  current  distribution  on  the  probe,  the  radial  current 
distribution  on  the  patch  and  the  azimuthal  magnetic  current  sheet  across  the  aperture  of 
the  driving  coaxial  line.  This  set  of  equations  is  then  solved  using  the  method  of  moments. 
The  resulting  matrix  equation  is  obtained  in  terms  of  Sommerfeld-type  integrals  that  take 
into  account  the  effect  of  the  layered  medium.  These  integrals  are  efficiently  computed  by 
a  simple  deformation  in  the  complex  wavenumber  domain.  The  probe  current  distribu¬ 
tion,  input  impedance  and  radiation  pattern  are  presented  and  compared  to  the  case  of  a 
uniform  probe  current  distribution. 

Microstrip  antennas  of  stacked  configurations  have  received  attention  in  recent  years 
for  both  wideband  and  dual  frequency  use,  overcoming  the  narrow  bandwidth  of  conven¬ 
tional  single  layer  microstrip  antennas.  Although  much  experimental  work  has  been  per¬ 
formed,  theoretical  analyses  of  stacked  microstrip  patches  is  limited.  Resonant  frequencies 
of  the  stacked  microstrip  antennas  have  been  rigorously  calculated.  Numerical  methods 
have  been  used  to  calculate  the  current  and  radiation  fields  of  a  stacked  microstrip  antenna. 
The  method  of  moments  has  been  applied  to  analyze  to  the  stacked  microstrip  structure 
when  excited  by  an  incident  plane  wave.  A  spectral  domain  iterative  analysis  for  a  stacked 
microstrip  antenna  where  the  antenna  is  described  by  a  rectangular  sampling  grid  has  been 
used  to  calculate  radiation  patterns.  This  analysis  does  not  allow  for  accurate  modeling 
of  the  probe  feed. 

In  our  approach,  the  input  impedance  and  radiation  fields  of  a  coaxial  probe-fed 
microstrip  antenna  consisting  of  two  circular  microstrip  disks  in  a  stacked  configuration 
are  investigated.  Using  a  dyadic  Green’s  formulation,  a  rigorous  analysis  of  the  microstrip 
antenna  is  performed  for  two  stacked  configurations.  Assuming  uniform  current  along  the 
probe,  the  mixed  boundary  value  problem  is  reduced  to  a  set  of  coupled  vector  integral 
equations  using  the  vector  Hankel  transform  and  solved  using  Galerkin’s  method  in  the 
spectral  domain.  Due  to  the  singular  nature  of  the  current  on  the  driven  disk  in  the  vicinity 
of  the  probe,  an  additional  term  is  included  in  the  current  expansion  to  account  for  the 
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divergent  nature  of  the  current  near  the  probe  feed  junction  and  insure  continuity  of  the 
current  at  the  junction. 

The  input  impedance  and  radiation  patterns  of  the  stacked  microstrip  antenna  is 
calculated  as  a  function  of  the  layered  substrate,  permittivities  and  thicknesses,  and  the 
ratio  of  the  radii  of  the  two  disks.  Both  dual  frequency  and  wideband  operation  is  discussed. 

Microstrip  discontinuities,  such  as  open  end,  gap  and  step  in  width,  have  been 
widely  studied  by  many  authors.  There  axe  different  methods  for  analyzing  microstrip  dis¬ 
continuities,  such  as  quasi-static  approach,  planar  waveguide  model  and  integral  equation 
formulation.  As  the  frequrncy  gets  higher,  the  quasi-static  assumption  is  not  valid.  In  the 
planar  waveguide  model  analyses,  the  thickness  of  the  substrate  is  assumed  much  smaller 
than  the  wavelength  so  that  a  two-dimensional  model  may  be  applied.  In  this  case,  the 
effect  of  the  radiation  and  the  surface  waves  axe  not  considered.  The  integral  equation 
method  has  been  applied  to  study  the  open  end  and  gap  dicontinuities  on  isotropic  sub¬ 
strates.  In  applying  the  integral  equation  method,  various  approximation  were  introduced 
in  the  computation  procedure.  More  recently,  finite  element  expansion  currents  are  used 
to  formulate  a  full- wave  analysis  of  micristrip  discontinuities  on  isotropic  substrates. 

The  open  end,  gap  and  step  in  width  discontinuities  placed  on  anisotropic  substrates 
are  rigorously  analyzed.  Both  uniaxial  and  tilted  uniaxial  anisotropy  are  considered.  The 
materials  are  assumed  to  be  lossless  and  the  metal  strips  to  be  infinitely  thin.  A  dyadic 
Green’s  function  for  layered  anisotropic  media  is  used  to  formulate  a  set  of  vector  integral 
equations  for  the  current  distribution.  The  fundamental  hybrid  mode  is  assumed  to  be 
propagating  on  the  input  and  output  of  microstrip  lines.  In  solving  the  set  of  vector  integral 
equations,  the  method  of  moment  is  employed.  The  basis  functions  for  the  current  on  the 
metal  strip  consider  the  edge  effect.  Both  logitudinal  said  transverse  currents  are  considered 
in  the  calculation.  The  propagation  constant  for  the  infinitely  long  uniform  microstrip  line 
is  first  calculated.  Then  the  propagation  constant  of  the  fundamental  mode  is  used  to 
formulate  the  excitation  of  the  discontinuity  problem.  At  the  discontinuity,  local  basis 
functions  are  used  to  simulate  the  local  currents  near  the  discontinuity.  The  scattering 
matrix  can  then  be  obtained,  and  an  equivalent  circuit  model  can  be  proposed.  The  effect 
of  the  anisotropy  is  investigated  and  the  results  are  discussed. 
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The  leakage  phenomenon  is  important  in  the  area  of  millimeter- wave  integrated 
circuits  and  integrated  optics.  Theoretical  analyses  and  experiments  have  been  performed 
to  investigate  this  phenomenon.  The  leakage  is  due  to  the  TE-TM  coupling  occurring  at 
the  geometrical  discontinuities,  and  the  leaky  power  in  the  form  of  surface  wave  propagates 
in  the  background  medium. 

There  are  different  methods  to  analyze  the  dielectric  strip  waveguides,  including 
the  approximate  field  matching  method,  effective  dielectric  constant  (EDC)  method,  mode 
matching  method,  etc.  The  first  two  methods  axe  approximate,  and  can  not  be  used  to 
predict  the  imaginary  part  of  the  propagation  constant.  In  the  third  one,  ground  planes 
have  to  be  put  at  some  distance  away  from  the  guiding  structure,  hence  the  effect  of 
radiation  loss  is  neglected. 

An  integral  equation  formulation  using  dyadic  Green’s  function  is  derived  to  solve 
for  the  dispersion  relation  of  single  and  coupled  dielectric  strip  waveguides.  A  method  to 
predict  the  leakage  is  presented,  and  the  leakage  properties  are  investigated. 

Three  different  dielectric  strip  waveguides  are  investigated  :  optical  rib  waveguide, 
strip  dielectric  guide,  and  insulated  image  guide.  Both  single  and  coupled  strip  waveguides 
are  studied.  The  cross  section  of  the  dielectric  strips  are  assumed  to  have  rectangular  shape. 
Applying  the  Galerkin’s  method,  the  field  distribution  on  the  cross  section  are  represented 
by  a  set  of  unit  pulse  basis  functions.  Substituting  these  basis  functions  into  the  integral 
equations,  and  choosing  the  same  set  of  basis  functions  as  the  testing  functions,  we  can 
obtain  a  determinant  equation  from  which  the  propagation  constant  can  be  solved. 

For  single  dielectric  strip  waveguide,  it  is  observed  that  the  leakage  occurs  when  the 
effective  refractive  index  is  smaller  than  that  of  a  surface  wave  mode  in  the  background 
medium.  It  is  also  observed  that  if  the  lowest  TE-like  (TM-like)  mode  is  leaky,  the  lowest 
TM-like  (TE-like)  mode  is  non-leaky.  When  the  lowest  order  mode  leaks,  the  surface  wave 
mode  of  opposite  polarization  is  excited.  When  the  higher  order  mode  leaks,  the  surface 
wave  modes  of  both  polarizations  can  be  excited. 


12 


For  two  symmetrical  dielectric  strip  waveguides,  both  the  even  and  odd  modes 
are  investigated.  For  the  leaky  mode,  the  total  leakage  is  due  to  the  leakage  from  each 
individual  strip  waveguide.  At  the  separation  where  the  even  mode  has  a  maximum  leakage, 
it  implies  that  the  surface  wave  modes  excited  by  each  waveguide  add  in  phase.  For  the 
odd  mode  at  about  the  same  separation,  these  coaxial  line  feed,  the  reflection  coefficient 
for  the  TEM  mode  is  obtained  which  allows  one  to  compute  the  input  impedance  at  the 
terminals  of  the  probe.  Numerical  results  for  the  input  impedance  are  presented. 

A  finite  difference  time  domain  technique  for  two  dimensional  time  domain  scat¬ 
tering  of  electromagnetic  waves  is  derived.  The  triangular  grids  and  the  control  region 
approximation  are  employed  to  discretize  Maxwell’s  equations.  The  finite  difference  time 
domain  techniques  with  uniform  rectangular  grids  has  been  used  in  the  past.  The  scatter¬ 
ed  are  modeled  using  staircases  and,  recently,  the  accuracy  of  this  approximation  has  been 
investigated.  Several  types  of  other  grids  have  been  proposed  to  improve  the  staircase  ap¬ 
proximation.  Generalized  nonorthogonal  grid  can  model  scatterer  without  staircasing.  It 
has  been  applied  to  spherical  systems,  yet  they  appear  to  be  cumbersome  for  general  scat¬ 
tered.  The  “distorted  rectangular  grid”  model  approximates  the  computational  domain 
using  rectangular  grids  and  distorts  the  boundary  grids  to  fit  the  interfaces.  The  triangular 
grid  is  used  in  this  paper,  which  is  very  flexible  in  dealing  with  arbitrary  scattered  and 
absorbing  boundaries. 

The  control  region  approximation,  which  calls  for  Delaunay  and  Dirichlet  tessella¬ 
tion,  has  been  successfully  applied  to  the  frequency  domain  problems  in  the  past.  Two 
double  integral  terms  are  obtained  by  integrating  the  Helmholtz  equation  about  the  De¬ 
launay  tessellation.  The  term  involving  the  Laplace  operator  can  be  converted  to  a  closed 
loop  integral  of  normal  derivatives,  which  can  easily  be  approximated  in  finite  difference 
manner  by  utilizing  the  orthogonal  property  of  Delaunay  and  Dirichlet  tessellation.  The 
remaining  term  can  be  approximated  by  multiplying  the  field  at  the  node  with  the  area.  In 
the  time  domain  problem,  the  same  approximation  is  applied  to  the  wave  equation,  except 
the  term  involving  time  derivatives  is  used  in  time  marching  scheme.  Alternatively,  as  in 
Yee’s  algorithm,  the  first  order  Maxwell’s  equations  are  solved  by  spatially  and  temporally 
separating  the  electric  and  magnetic  fields.  In  the  case  of  electric  polarization,  the  electric 
fields  are  placed  at  the  nodes  and  the  magnetic  fields  are  placed  at  the  center  of  triangular 
edges.  The  curl  H  equation  is  integrated  by  applying  Stoke’s  theorem  and  convert  it  to  a 
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closed  loop  integral  of  tangential  magnetic  fields.  This  equation  can  be  used  to  advance 
electric  fields  in  time.  To  update  magnetic  fields,  the  second  curl  equation  is  used.  This 
equation  is  approximated  in  the  finite  difference  manner  by  utilizing  the  orthogonality 
property  of  the  tessellation.  The  equations  for  the  magnetic  polarization  case  can  also  be 
derived  following  the  similar  procedure. 

In  order  to  limit  the  computation  domain,  the  scatterers  are  enclosed  with  artificial 
outer  boundaries.  Continuous  smooth  outer  boundaries,  such  as  circles  and  ellipses,  are 
chosen.  The  second-order  time  domain  absorbing  boundary  conditions  derived  from  the 
pseudo-differential  operator  approach  is  imposed  at  the  outer  boundaries.  These  boundary 
conditions  are  implemented  with  the  control  region  approximation  to  determine  necessary 
field  quantities  at  the  boundary.  The  results  of  the  time  domain  control  region  approach 
are  presented  for  simple  scatterer  geometries,  such  as  conducting  and  coated  cylinders  and 
strips,  by  calculating  both  the  transient  and  time-harmonic  responses. 

The  Finite-Difference  Time-Domain  (FD-TD)  method  was  first  introduced  by  Yee 
who  discretized  Maxwell’s  time  dependent  curl  equations  with  second-order  accurate 
central-difference  approximations  in  both  the  space  and  time  derivatives.  Since  then,  it 
has  been  applied  extensively  to  scattering  and  wave  absorption  problems.  Application  of 
the  FD-TD  method  to  microstrip  problems,  in  which  frequency-domain  approaches  have 
dominated,  has  so  far  attracted  little  attention  until  recently  it  was  used  to  obtain  fre¬ 
quency  characteristics  of  microstrip  cavities.  Also,  it  has  been  extended  to  the  analyses 
of  open  microstrip  line  and  microstrip  discontinuity  problems  where  absorbing  boundary 
conditions  are  needed  for  the  simulation  of  the  unbounded  domain.  However,  only  isotropic 
or  simple  anisotropic  media  are  considered  in  the  above  papers. 

A  new  FD-TD  grid  model  is  used  to  solve  microstrip  problems  in  anisotropic  media 
having  tilted  optical  axes  expressed  by  permittivity  or  permeability  tensor  with  off-diagonal 
elements.  This  grid  model  is  indeed  a  superposition  of  two  conventional  grids  with  some 
displacement  which  depends  on  the  optical  axes  of  anisotropy.  Implementations  of  different 
boundary  conditions  are  discussed.  Using  this  model,  the  frequency-dependent  character¬ 
istics  of  microstrip  lines  are  investigated.  The  microstrips  are  assumed  to  be  placed  on 
top  of  anisotropic  substrates  with  tilted  optical  axes.  The  case  with  superstrates  is  also 
investigated. 
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In  the  finite  difference  computation,  the  open-end  termination  is  simulated  by  using 
the  open-circuit,  short-circuit  technique.  The  source  plane  is  implemented  by  using  a 
magnetic  wall  with  a  Gaussian  pulse  excited  on  the  surface  under  the  strip.  Because  of  the 
symmetry  of  the  problem,  the  region  under  consideration  can  be  reduced  by  half,  using  a 
magnetic-wall  at  the  center  plane. 

The  fields  at  different  positions  are  first  calulated.  Then  the  Fourier  Transform  is 
taken  to  give  the  field  spectra  from  which  the  voltage  find  current  can  also  be  obtained. 
Using  these  data,  the  effective  permittivity  and  the  characteristic  impedance  can  be  de¬ 
termined.  The  frequency  characteristics  of  microstrip  lines  in  anisotropic  media  obtained 
by  this  method  are  compared  with  the  published  results. 

Finite  difference  time  domain  (FDTD)  techniques  show  great  promise  in  their  ability 
to  solve  three  dimensional  problems  with  arbitrary  geometry.  Advantages  of  this  method 
include  the  ability  to  model  spatially  or  temporally  varying  media.  These  advantages 
are  due  to  the  complete  discretization  of  both  space  and  time.  Considering  the  volume 
of  information  being  calculated  these  techniques  are  very  efficient  and  are  well  suited  to 
calculation  on  future  parallel  processing  computers.  This  method  was  first  formulated 
by  Yee  in  1966  and  his  basic  algorithm  is  still  in  use.  Recent  work  has  demonstrated  the 
applicability  of  the  FDTD  technique  to  microstrip  problems.  The  centered  finite  difference 
approximations  used  are  second  order  accurate  in  both  space  and  time  yielding  good  results 
for  reasonable  mesh  sizes.  Numerical  techniques  used  to  solve  electromagnetic  problems 
must  limit  the  domain  over  which  the  fields  are  to  be  calculated.  This  mandates  the  use  of 
an  absorbing  boundary  condition  to  simulate  the  outward  propagation  of  waves  incident 
on  the  walls  of  the  mesh.  An  absorbing  boundary  condition  has  been  developed  by  Mur 
based  on  the  work  of  Enquist  and  Majda. 
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Our  work  in  this  area  includes  development  of  the  algorithms  mentioned  above  into 
a  general  purpose  computer  code  which  may  be  used  to  solve  for  the  transient  response 
of  electromagnetic  problems  with  an  arbitrary  geometry.  In  addition  to  the  transient 
response,  frequency  domain  parameters  may  be  obtained  by  fourier  transform  of  the  time 
domain  results.  Since  the  fields  are  calculated  throughout  space  and  time  all  other  desired 
parameters  may  be  calculated  from  the  field  quantities.  Specifically,  we  are  analyzing 
rectangular  microstrip  structures  with  as  many  as  two  or  more  ports.  Such  structures  may 
be  used  in  MMIC  filters  or  antennas.  This  problem  is  of  interest  for  several  reasons.  First, 
there  are  existing  frequency  domain  solutions  to  the  resonance  problem  of  a  rectangular 
microstrip  patch,  which  we  may  compare  with  the  FDTD  solution.  Secondly,  the  FDTD 
technique  may  be  used  to  analyze  coupling  of  microstrip  lines  to  the  rectangular  structure. 
This  coupling  may  be  either  a  direct  connection  or  a  gap  coupled  connection.  Advantages 
of  the  FDTD  solution  of  this  problem  are  that  it  is  a  full  wave  solution  which  allows  for 
radiation  or  surface  wave  loss  and  that  no  empirical  values  such  as  "effective”  dimensions 
are  needed  for  the  analysis,  also  the  geometry  may  be  altered  easily  to  allow  for  various 
connections  or  coupling  to  the  patch.  This  is  a  significant  improvement  over  methods 
which  rely  on  a  planar  circuit  approach  in  which  the  substrate  thickness  must  be  small 
compared  to  wavelength  and  inherently  three  dimensional  coupling  problems  are  not  easily 
handled.  Comparison  of  our  results  with  various  planar  circuit  approaches  will  be  made. 

A  new  perturbation  series,  coupled  integral  equation  approach  for  calculating  the 
frequency  dependent  circuit  parameters  for  quasi-TEM  transmission  lines  with  lossy  con¬ 
ductors  is  presented.  The  method  considers  the  addition  of  loss  and  dispersion  to  be 
perturbations  on  the  lossless  TEM  case,  and  therefore  the  difference  between  the  propaga¬ 
tion  constant  and  the  wavenumber  in  free  space  is  a  small  parameter.  We  obtain  the  lowest 
order  term  of  the  perturbation  series  by  solving  two  quasistatic  problems;  the  electrostatic 
problem  to  get  the  capacitance,  and  the  magnetoquasistatic  problem,  with  the  distribution 
of  current  inside  the  wire  considered,  which  gives  the  frequency-dependent  inductance  and 
resistance.  Both  of  these  problems  are  solved  using  one-dimensional  integral  equations 
for  quantities  on  the  surface  of  the  conductor;  this  represents  a  significant  improvement 
in  efficiency  over  previous  methods.  For  most  cases  of  practical  interest,  the  lowest  order 
term  of  the  series  will  suffice.  If,  however,  the  change  in  the  propagation  constant  from  the 
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lossless  case,  due  to  the  altered  inductance  and  the  addition  of  resistance,  is  significant, 
additional  terms  in  the  perturbation  series  can  be  calculated. 

The  method  is  illustrated  with  the  case  of  one  or  more  wires  embedded  in  a  uniform 
dielectric.  In  the  original  magnetoquasistatic  problem,  the  current  is  entirely  directed  along 
the  axis  of  propagation,  and  satisfies  the  frequency-domain  diffusion  equation.  Outside  the 
wire,  the  magnetic  vector  potential  is  in  the  same  direction,  and  obeys  Laplace’s  equation. 
The  boundary  conditions  are  the  continuity  of  tangential  and  normal  magnetic  field  at  the 
interface,  which  can  be  expressed  in  terms  of  the  current  density  and  vector  potential  and 
their  derivatives.  Since  we  can  express  the  ratio  of  the  frequency-dependent  resistance  to 
the  DC  resistance  in  terms  of  the  values  of  the  volume  current  and  its  normal  derivative 
on  the  surface  of  the  wire  only,  we  can  use  a  pair  of  coupled  integral  equations  to  solve  for 
these  quantities  alone,  which  we  can  solve  by  Galerkin’s  method  or  other  finite  element 
methods. 

Results  obtained  using  this  technique  are  shown  for  some  important  cases,  includ¬ 
ing  rectangular  wires,  and  are  compared  with  earlier  methods  and  with  experimental  data. 
Previous  methods  for  calculating  the  resistance  fall  into  three  categories.  First,  for  cer¬ 
tain  cases,  exact  analytical  results  can  be  obtained.  Secondly,  especially  in  the  case  of  a 
rectangular  wire,  the  cross-section  can  be  divided  into  rectangular  segments,  each  much 
smaller  than  a  skin-depth,  across  which  the  current  is  assumed  to  be  constant.  Magne¬ 
toquasistatics  gives  simple  answers  for  the  resistance  and  self-inductance  of  each  element, 
and  the  mutual  inductance  between  elements.  This  leads  to  a  matrix  equation  which  is 
solved  for  the  current  distribution.  The  disadvantages  of  this  technique  are  that  it  requires 
basis  functions  throughout  the  cross-section  of  the  conductor,  which  is  especially  intensive 
as  the  frequency  gets  large.  Also,  closed  form  expressions  for  the  matrix  elements  only 
exist  when  the  elements  are  rectangular  —  other  shapes,  such  as  triangular  patches,  which 
might  be  used  to  fit  a  wire  of  arbitrary  shape,  lead  to  nested  numerical  integrals. 
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The  third  method  used  is  a  variational  procedure.  This  is  similar  to  the  method 
presented  here,  except  that  the  current  and  the  magnetic  vector  potential  are  expanded  in 
functions  which  span  the  entire  cross-section.  This  has  two  drawbacks:  first,  it  requires  that 
there  be  a  closed  outer  conductor,  which  is  not  physical  in  many  important  cases.  Second, 
as  in  the  previous  method,  using  elements  which  fill  the  entire  cross-section  increases 
the  computation  time  unnecessarily,  since  only  the  value  of  the  current  and  its  normal 
derivative  at  the  surface  of  the  wire  are  needed  to  calculate  the  resistance. 

A  new  method  for  analyzing  frequency-dependent  transmission  line  systems  with 
nonlinear  terminations  is  presented.  The  generalized  scattering  matrix  formulation  is  used 
as  the  foundation  for  the  time  domain  iteration  scheme.  Compared  to  the  admittance 
matrix  approach  proposed  in  a  previous  paper,  it  has  the  advantage  of  shorter  impulse 
response  which  leads  to  smaller  computer  memory  requirement  and  faster  computation 
time.  Examples  of  a  microstrip  line  loaded  with  nonlinear  elements  are  given  to  illustrate 
the  efficiency  of  this  method. 

As  the  speeds  of  integrated  circuits  increase,  the  effect  of  interconnection  lines  be¬ 
comes  more  and  more  important.  Traditional  lumped  element  circuit  model  must  be 
supplemented  by  the  transmission  line  model  in  order  to  account  for  propagation  de¬ 
lays,  dispersion  and  losses.  This  has  created  needs  for  new  numerical  procedures  that 
can  be  easily  incorporated  into  current  CAD  tools.  To  make  matters  more  complicated, 
the  interconnection  lines  are  terminated  with  not  only  linear  resistors  but  also  nonlinear 
semiconductor  devices,  such  as  diodes  and  transistors. 

Several  techniques  are  now  commonly  used  to  deal  with  nonlinear  circuit  prob¬ 
lems,  for  example,  the  direct  time  domain  approaches,  and  the  semi-frequency-domain  ap¬ 
proaches,  such  as  the  harmonic  balance  technique.  The  semi-frequency-domain  approaches 
are  useful  for  microwave  and  millimeter  wave  integrated  circuits  but  become  impractical 
for  digital  integrated  circuits  because  of  the  latter’s  wide  band  nature.  On  the  other  hand, 
frequency-dependent  parameters  often  make  it  awkward  to  apply  the  direct  time  domain 
approach  to  the  interconnection  line  systems. 
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We  propose  a  hybrid  frequency-domain  time-domain  technique  based  on  the  gen¬ 
eralized  scattering  matrix  formulation.  For  an  n  line  system,  we  define  2 n  scattering 
parameters  according  to  the  frequency-dependent  characteristic  impedances  of  individual 
lines  (Zqj  =  ^JLjj(u)/Cjj(u; )).  The  time-domain  transfer  matrix  (impulse  response)  of 
this  2n-port  system  is  then  obtained  by  the  inverse  Fourier  Transform.  Lastly,  the  non¬ 
linear  equations  associated  with  terminal  characteristics  are  incorporated  and  solved  with 
iteration  procedures  such  as  the  Newton- Ralphson  method. 

The  key  to  efficient  and  stable  solutions  in  this  problem  is  shortening  the  duration 
of  every  transfer  matrix  element.  With  the  generalized  scattering  parameters  approach,  we 
are  able  to  achieve  that  yet  eliminating  the  need  for  artificial  matching  networks  adopted 
by  a  previous  work.  Furthermore,  the  use  of  individual  characteristic  impedances  in  the 
definition  of  scattering  parameters  enables  us  to  generalize  this  method  to  coupled  lines 
with  distinct  properties  while  keeping  the  duration  of  transfer  matrix  elements  short.  This 
cannot  be  realized  if  traditional  scattering  parameters  are  used.  We  shall  illustrate  the 
elegance  and  efficiency  of  our  approach  for  a  dispersive  microstrip  line  with  different  non¬ 
linear  loads  and  excited  with  narrow  Gaussian  pulses.  The  elements  of  transfer  matrix 
are  found  out  to  be  either  zero  or  single  retarded  delta-impulse  accompanying  a  small 
spike  with  very  narrow  spread.  Typical  computation  time  for  a  1000  time-step  iteration 
ranges  from  4  to  27  seconds  on  a  VAXStation  3500.  The  effects  of  line  dispersion  and  load 
nonlinearity  will  be  clearly  delineated  in  the  presentation. 

The  transient  propagation  characteristics  of  VLSI  interconnects  with  discrete  ca¬ 
pacitive  loads  at  various  locations  is  analyzed  based  on  a  hybrid  transmission  lines-lumped 
element  circuit  model.  Exact  expressions  of  the  Laplace  transform  of  unit  step  responses 
are  first  obtained  through  the  ABCD  matrix  formulation.  We  then  apply  the  equivalent 
dominant  pole  approximation  to  the  transfer  function  with  the  propagation  delays  factored 
out.  The  approximated  transfer  function  can  be  inverted  in  closed  form  and  quickly  eval¬ 
uated.  These  results  provide  efficient  ways  of  finding  approximately  the  effects  on  delays 
and  rise  time  brought  by  VLSI  off-chip  interconnects. 
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Because  of  the  dramatic  increase  in  device  densities  on  microelectronic  chips,  the 
propagation  delay  for  off-chip  interconnects  has  become  the  limiting  factor  to  the  speed 
of  VLSI  packages.  Typical  scales  of  these  interconnects  will  be  comparable  or  larger  to 
the  characteristic  wavelength  of  high  frequency  components  of  the  signal.  Therefore,  to 
calculate  the  delays  caused  by  these  interconnects  properly,  a  hybrid  circuit  model  con¬ 
taining  transmission  line  sections  as  well  as  lumped  elements  must  be  used  in  place  of 
the  all-lumped  element  one.  Most  circuit  simulation  packages  are  nevertheless  based  on 
the  latter  and  have  to  resort  to  subsection  approximation  when  dealing  with  transmission 
lines.  This  scheme  will  undoubtedly  lead  to  lengthy  computation,  which  is  not  desirable 
when  a  quick,  heuristic  estimate  of  bounds  are  needed  for  the  initial  phase  of  the  design 
cycle. 

Two  approaches  have  been  developed  for  obtaining  the  approximate  transient  re¬ 
sponse  without  lengthy  simulation.  The  first  kind  of  solution  techniques  emphasize  the 
calculation  of  bounds  to  voltage  responses  from  the  differential  equations  either  by  direct 
integration  or  by  using  the  optimal  control  theory.  On  the  other  hand,  the  second  kind  of 
techniques  analyze  the  properties  of  Laplace  transform  domain  solution.  Thus  far,  their  ap¬ 
plications  are  limited  to  all  lumped-element  and  distributed  RC  networks,  which  can  only 
take  care  of  on-chip  interconnects.  We  shall  take  the  second  approach  by  incorporating 
transmission  line  elements  for  off-chip  delay  estimation. 

Our  configuration  includes  a  series  of  transmission  line  sections  with  arbitrary  dis¬ 
crete  capacitances  and  resistances  loaded  at  junctions.  The  ABCD  matrix  formulation  is 
used  to  obtain  the  Laplace  transform  of  the  unit  step  response.  We  express  the  latter  in  the 
form  of  exp(— sT)/Q(s),  where  Q(s)  =  Aq(j)  +  Ai(s)  exp(sri)  +  ^2(3)  exp(sT2)  +  •  •  •  with 
all  {Aj(s)}  being  rational  functions  in  s.  The  factor  {exp(— sT)}  is  identified  with  direct 
transmission  delay  over  the  total  length  of  the  line.  For  the  rest  part  (l/(J(s)),  we  proceed 
to  apply  the  equivalent  dominant  pole  approximation  technique[7].  Either  a  single  negative 
real  pole  or  second-order  complex  conjugate  pair  will  be  chosen  for  approximation  depend¬ 
ing  upon  the  property  of  lumped  loads  at  junctions  as  well  as  the  source  impedances.  A 
phase-correction  factor  exp(— aTm)  is  introduced  to  make  up  for  the  discrepancies  caused 
by  our  low-order  approximation.  The  first-order  and  second  order  approximations  enable 
us  to  obtain  closed-form  solution  to  the  transient  response.  Comparison  of  the  approxi¬ 
mated  responses  with  those  obtained  from  brute-force  numerical  Laplace  inversion  shows 
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very  good  match  when  the  propagation  delay  of  an  average  transmission  line  section  is 
less  than  half  the  product  of  junction  load  capacitance  and  transmission  line  characteristic 
impedance.  Yet  we  only  have  to  spend  a  fraction  of  the  time  for  computations.  The  accu¬ 
racy  of  this  method  will  be  discussed  in  detail  with  some  examples  of  lossless  transmission 
line  networks  in  which  lumped  capacitors  are  loaded  at  regular  intervals. 

The  propagation  properties  of  single  and  coupled  inhomogeneous  slab  waveguides 
are  analyzed.  An  integral  equation  formulation  using  the  dyadic  Green’s  function  which 
covers  both  the  TE  and  TM  modes  is  proposed.  The  dispersion  relations  are  obtained  by 
applying  the  Galerkin’s  method  to  solve  the  integral  equation.  The  coupling  between  two 
symmetrical  inhomogeneous  slab  waveguides  is  also  investigated.  This  method  is  shown 
to  be  applicable  to  arbitrary  dielectric  constant  profiles. 

The  guidance  and  leakage  properties  of  single  and  coupled  dielectric  strip  waveg¬ 
uides  are  analyzed  using  the  dyadic  Green’s  function  and  integral  equation  formulation. 
Galerkin’s  method  is  used  to  solve  the  integral  equation  for  the  dispersion  relation.  The 
effects  of  the  geometrical  and  the  electrical  parameters  on  the  dispersion  relation  are  inves¬ 
tigated.  A  method  to  predict  the  occurrence  of  leakage  is  proposed.  The  properties  of  the 
even  and  the  odd  leaky  modes  are  also  investigated.  Results  are  compared  with  previous 
analysis  and  shown  to  be  in  good  agreement. 

A  spectral  domain  dyadic  Green’s  function  for  multilayered  uniaxially  anisotropic 
media  containing  three-dimensional  sources  is  derived.  Tractable  forms  are  shown  to  be 
easily  deduced  from  the  physical  picture  of  the  waves  radiated  by  the  primary  sources 
and  the  multiple  reflections  from  the  stratified  medium.  The  formulation  decomposes 
the  dyadic  Green’s  function  into  TE  and  TM  waves.  The  dyadic  Green’s  function  in  the 
source  region  is  properly  represented  by  extracting  the  delta  function  singularity.  A  simple 
proceedure  to  obtain  the  fields  in  any  arbitrary  layer  is  described.  Recursion  relations  for 
appropriately  defined  reflection  and  transmission  coefficients  are  presented.  Forms  suitable 
for  transmission  line  applications  in  multilayered  media  are  derived. 
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Full  modal  analysis  is  used  to  study  the  dispersion  characteristics  of  microstrip 
lines  periodically  loaded  with  crossing  strips  in  a  stratified  uni  axially  anisotropic  medium. 
Dyadic  Green’s  functions  in  the  spectral  domain  for  the  multilayered  medium  in  conjunc¬ 
tion  with  the  vector  Fourier  transform  (VFT)  are  used  to  formulate  a  coupled  set  of  vector 
integral  equations  for  the  current  distribution  on  the  signal  line  and  the  crossing  strips. 
Galerkin’s  procedure  is  applied  to  derive  the  eigenvalue  equation  for  the  propagation  con¬ 
stant.  The  effect  of  anisotropy  for  both  open  and  shielded  structures  on  the  stopband 
properties  is  investigated. 

A  direct  three  dimensional  finite  difference  time  domain  (FDTD)  method  is  applied 
to  the  full-wave  analysis  of  various  microstrip  structures.  The  method  is  shown  to  be  an 
efficient  tool  for  modelling  complicated  microstrip  circuit  components  as  well  as  microstrip 
antennas.  From  the  time  domain  results,  the  input  impedance  of  a  line-fed  rectangular 
patch  antenna  and  the  frequency  dependent  scattering  parameters  of  a  low  pass  filter  and 
a  branch  line  coupler  are  calculated.  These  circuits  are  fabricated  and  the  measurements 
are  compared  with  the  FDTD  results  and  shown  to  be  in  good  agreement. 

A  rigorous  dyadic  Green’s  function  formulation  in  the  spectral  domain  is  used  to 
study  the  dispersion  characteristics  of  signal  strip  lines  in  the  presence  of  metallic  crossing 
strips.  A  set  of  coupled  vector  integral  equations  for  the  current  distribution  on  the 
conductors  is  derived.  Galerkin’s  method  is  then  applied  to  derive  the  matrix  eigenvalue 
equation  for  the  propagation  constant.  The  dispersion  properties  of  the  signal  lines  are 
studied  for  both  cases  of  finite  and  infinite  length  crossing  strips.  The  effects  of  the 
structure  dimensions  on  the  passband  and  stopband  characteristics  are  investigated.  For 
crossing  strips  of  finite  length,  the  stopband  is  mainly  affected  by  the  period,  the  crossing 
strip  length,  and  the  separation  between  the  signal  and  the  crossing  strips.  For  crossing 
strips  of  infinite  length  carrying  travelling  waves,  attenuation  along  the  signal  line  exists 
over  the  whole  frequency  range  of  operation. 
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Abstract 

A  macroscopic  model  is  proposed  to  explain  nonlinear  electromagnetic  phenomena  in  su¬ 
perconductors.  Nonlinear  constitutive  relations  for  electromagnetic  problems  are  derived 
by  modifying  the  linear  London’s  equations.  Tne  superelectron  number  density  n3  is  a 
function  of  the  applied  current  density  J.  The  critical  current  density  Jc  is  derived  clas¬ 
sically  from  a  critical  energy  Ec.  For  temperature  2  ^  0,  the  concept  of  critical  current 
Jc  does  not  imply  an  abrupt  transition  of  the  whole  sample  from  a  superconducting  state 
to  a  normal  state  when  J  >  Jc.  A  rather  smooth  variation  of  ns(J)  is  shown  instead. 
The  relation,  na(  J),  is  derived  from  the  Maxwellian  distribution  of  electron  velocities  at 
a  certain  temperature  T  and  a  certain  macroscopic  current  density  J.  Agreement  has 
also  been  found  between  this  nj(/,T)  model  and  the  temperature  dependence  of  n$  in 
the  two-fluid  model.  The  nonlinear  conductivities  <rs(J)  and  <rn(J)  are  obtained  from  the 
London’s  equation  and  the  ns(J)  function.  Nonlinear  resistance  i2(J),  kinetic  inductance 
Lf.(I)  and  surface  impedance  Z3(I)  in  thin  wire,  slab,  and  strip  geometries  of  supercon- 
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ductors  axe  calculated.  A  general  scheme  of  solving  nonlinear  electromagnetic  problems  in 
superconductors  is  proposed.  A  good  agreement  between  the  theory  and  experiments  has 
been  found. 
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I.  Introduction 


Superconductors  have  great  potential  applications  in  many  fields.  For  example,  in  mi¬ 
crowave  integrated  circuits,  high-Q  resonators  and  microstrips  can  be  made  of  supercon¬ 
ductors  for  lower  losses.  The  discovery  of  high-Tc  superconductors  has  also  made  a  big 

G 

impact  on  superconductor  modeling.  The  problem  of  modeling  nonlinear  electromagnetic 
properties  of  a  superconductor  is  of  practical  importance.  For  example,  in  application 
of  superconductors  to  high-Q  resonators,  fields  and  currents  are  very  large  at  resonance. 
Nonlinear  effects  are  inevitable. 

In  this  paper,  we  use  the  macroscopic  (classical)  theory  to  model  the  nonlinear  supercon¬ 
ductivity.  The  constitutive  relations  which  relate  electric  field  E  and  magnetic  field  B 
to  superconducting  current  density  Js  will  be  derived.  The  theory  is  based  on  the  two 
London  equations  and  the  two-fluid  model.  The  new  problem  at  hand  is  to  incorporate  the 
nonlinear  effects  into  the  constitutive  relations.  Application  of  the  constitutive  relations 
will  provide  new  methods  for  studying  nonlinear  effects  in  superconductors. 

There  are  different  types  of  nonlinearity  in  superconductors.  One  type  of  nonlinearity  is 
displayed  in  polycrystalline  superconductors  where  granular  currents  are  involved  [1,2]. 
The  granular  nonlinearity  occurs  when  the  current  I  is  above  Ic  (here  Ic  is  the  threshold 
current  of  the  grain  junctions).  It  is  similar  to  a  p-n  junction’s  exponential  I-V  relation. 
It  is  also  found  in  type  II  superconductors  that  the  vortex  motion  in  the  mixing  state 
(between  the  superconducting  and  the  normal  states)  can  cause  a  nonlinear  V  —  I  relation 
when  current  density  J  is  slightly  greater  than  a  critical  current  density  Jc  [3,  p.331]  [4, 
p.373].  The  granular  or  vortex  nonlinearity  will  not  be  discussed  further  in  this  paper. 

Another  type  of  nonlinearity  is  intrinsic  for  all  superconductors.  This  is  due  to  the  de- 
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pendence  of  the  superelectron  number  density  ns  on  the  applied  current  density  J.  This 
nonlinearity  is  more  general.  There  have  been  few  papers  in  literature  directly  addressing 
this  nonlinearity  problem.  In  this  paper  we  will  focus  on  this  intrinsic  nonlinearity. 

The  purpose  of  this  paper  is  to  provide  engineers  a  material-independent  macroscopic 
model  for  nonlinearities  in  superconductors.  The  model  will  be  based  on  macroscopic 
parameters.  The  model  is  intended  to  achieve  the  following:  l]  to  explain  the  nonlinear 
voltage-current  ( V  —  I)  relations  and  the  dependence  of  inductances  on  currents  (L  —  I 
relations).  This  requires  the  derivation  of  a  complex  conductivity  ct(  J)  model;  2]  to  explain 
the  concept  of  critical  current  Jc  with  no  abrupt  transition  in  the  sample’s  V  —  I  curve 
when  J  exceeds  Jc;  3]  to  obtain  the  temperature  dependence  of  ns(T)  which  agrees  with 
what  is  assumed  in  the  two-fluid  model;  4]  to  derive  a  general  scheme  which  can  solve  the 
nonlinear  problem  including  the  effects  of  geometry. 

II.  Critical  Energy  for  a  Single  Electron,  EC{T) 

The  starting  point  of  this  model  assumes  that  there  is  a  critical  value  of  energy,  EC)  which 
separates  the  superconducting  and  normal  states  of  a  single  electron.  When  an  electron 
has  an  energy  E  <  Ec,  it  is  in  the  superconducting  state,  it  is  paired  with  another  electron; 
when  E  >  Ec,  it  is  in  the  normal  state.  This  critical  value  Ec  is  a  function  of  temperature 
T.  The  higher  T  is,  the  lower  Ec  becomes.  At  T  =  Tc,  the  critical  temperature,  Ec  =  0. 
This  critical  energy  Ec  reminds  us  the  gap  parameter  A  in  the  BCS  microscopic  theory.  It 
may  correspond  to  the  de-pairing  energy  of  a  Cooper  pair.  But  here  the  critical  energy  Ec 
is  proposed  as  merely  an  assumption  without  justification.  The  origin  of  this  Ec  should 
be  studied  in  microscopic  theory.  We  will  leave  this  problem  and  use  Ec  as  a  starting 
assumption. 
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We  assume  the  temperature  dependence  of  Ec  as 

EC(T)  =  Z.52kBTc(l  -  £-)at  for  0  <  T  <  Tc  (1) 

J-c 

where  kg  is  the  Boltzmann  constant,  Tc  is  the  critical  temperature,  and  a  is  a  free  pa¬ 
rameter  which  may  depend  on  the  type  of  material,  a  >  0  is  required. 

For  each  electron,  the  energy  consists  of  kinetic  energy  and  potential  energy, 

1  9 

E  =  Ek  ■+■  Ep  =  -mv  +  Ep  (2) 

where  m  is  the  mass  of  an  electron  and  v  is  its  velocity.  In  quantum  theory,  v  is  the 
expectation  value  <  ip\v\ ip  >,  where  ip  is  the  wavefunction  of  the  electron  and  v  is  the 
velocity  operator.  The  electron  energy  E  can  be  expressed  as  a  function  of  temperature 
T,  applied  current  density  J,  magnetic  field  H,  and  field  frequency  /. 

E  =  E(T,JJ,H)  (3) 

Ec  is  the  origin  for  critical  values  of  Tc,  Jc,  fc  and  JJC.  It  is  useful  to  determine  the 
relations  of  these  critical  values  to  Ec. 


HI.  Velocity  Distribution  for  a  Group  of  Electrons 


First  let’s  consider  that  the  electrons  have  one-dimensional  velocities,  for  example,  in  a 
thin  (radius  o  <<  A  the  penetration  depth)  wire  within  which  the  current  flows  in  only  one 
direction.  The  electrons  in  the  sample  have  different  velocities  due  to  thermal  motion.  The 
average  velocity  of  all  electrons  is  nonzero  along  the  current  direction.  The  current  density 
J  is  a  macroscopic  quantity.  At  a  certain  point  r,  J  is  related  to  the  average  velocity  of 
the  electrons  in  a  small  volume  AV , 

EjVi 
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where  N  is  the  number  of  electrons  in  Ay  and  q  is  the  charge  of  an  electron.  Note  that 
not  every  electron  has  the  same  velocity.  The  velocities  t>t-  of  the  electrons  obey  a  certain 
distribution.  Here,  for  a  particular  example,  we  assume  that  the  number  of  electrons, 
6Nt  which  have  the  velocities  between  v  and  6v  obeys  the  one-dimensional  Maxwellian 
distribution.  If  we  define 


n(v)  = 


SN 

SvAV' 


(5) 


then  [10,  p.80]  [11,  p.70] 

<6> 

where  is  the  average  velocity  of  the  electrons,  m  is  the  mass  of  an  electron,  and  n* 
is  a  constant.  The  distribution  is  Gaussian.  Figure  1  shows  the  n(v)  function.  n(v)  is 
assumed  to  be  a  continuous  distribution.  The  state  of  an  electron  will  be  determined  from 
the  velocity  (energy)  in  this  distribution  n(t>). 


The  constant  n*  can  be  determined  from  the  normalization  condition, 

*+oo 


I. 


oo 


n(v)dv  —  ng 


(7) 


where  tiq  is  the  total  number  density  of  electrons  known  at  AV.  It  can  be  easily  shown 
that 

nt  ~  n0-  (8) 


We  can  relate  the  two  macroscopic  quantities,  J  and  vj±.  It  can  be  shown  that 

*+oo 


/+°0 

qn(v)vdv  =  qriQVj^ 

-00 


(9) 


at  Ay. 


The  T  dependence  of  n(v)  shows  that  at  T  =  0, 


n(v)  =  no$(«)> 
6 


(10) 


ft  6  function  [12,  p.16].  When  T  becomes  larger,  the  width  y/2kgTfm  of  the  distribution 


is  wider. 


IV.  ns(J)  and  Critical  Current  Density  Jc(T) 

Since  the  electrons  have  different  velocities,  therefore  different  energies  at  a  certain  applied 
current  density  J,  the  electrons  will  not  all  exceed  the  critical  energy  Ec  at  the  same  time. 
There  is  no  macroscopic  critical  value  Jc  of  the  current  density  such  that  when  J  >  Jc, 
number  density  of  superelectrons  na  =  0.  Only  when  T  =  0  and  n(v)  becomes  nQ$(v), 
na(  J)  shows  an  abrupt  drop  to  zero  at  a  critical  JC(T  =  0):  for  J  <  Jc(0),  na  =  ng;  and  for 
J  >  Jc( 0),  na  =  0.  In  experiments,  since  absolute  T  =  OK  is  not  achievable,  a  completely 
sharp  transition  of  na  is  not  observed. 

We  can  define  a  Jc(T)  which  corresponds  to  the  average  velocity  and  Ec.  First,  let’s 
define  a  critical  velocity  vc  for  a  single  electron  corresponding  to  Ec.  If  the  potential  energy 
is  included  in  the  case  of  maximum  kinetic  energy,  jmt>c  =  Ec-  Hence, 

”e = <n> 

Since  J  =  nogw^,  we  define 


Jc(T)  -  qn0vc  -  qn0J  J  ^ 

(12) 

MT)  =  -  y/'2 

(13) 

as  the  critical  current  density.  This  Jc  is  different  from  the  Ginzburg-Landau’s  de-pairing 
current  density.  Here  Jc(T)  is  a  quantity  derived  from  the  hypothetical  critical  energy  Ec 
via  the  classical  kinetic  energy  expression.  The  concept  of  this  critical  current  density  Jc 
is  also  different. 
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At  T  ^  0,  n3(  J)  is  a  smoothly  varying  function.  The  Jc{T)  has  a  different  characteristics 
from  Jc(T  =  0).  When  J  >  Jc{T),  n3  does  not  go  to  zero.  This  smooth  varying  feature  of 
the  I  —  V  curve  has  been  widely  observed  in  experiments  [5,6]. 

Since  =  J/qnQ,  the  number  density  of  superelectrons  (here  we  count  the  single  electron 
density  rather  than  the  pair  density)  na(J)  can  be  derived  from  the  Maxwellian  distribu¬ 
tion.  The  superelectrons  are  those  whose  energy  are  lower  than  Ec. 

mj) = l-i n(v)dv = C2)  n° («) 

In  general,  this  integral  can  only  be  evaluated  numerically.  A  special  case  is  at  T  =  0, 
n(v)  =  ng£(v  — 1>^).  Hence,  ns(J )  =  ng  if  <  vc  or  J  <  Jc( 0);  ns( J)  =  0  if  vj±  >  vc  or 
J  >  Jc( 0). 

It  is  illustrated  in  Fig.  2  that  when  Vji  +  y/2kj}T/m  <  vc,  almost  all  electrons  are 
superelectrons,  hence  n3  «  ng;  when  —  y/2 kjjT/m  >  vc,  almost  all  electrons  axe 
normal  electrons,  hence  ns  «  0.  At  T  =  Te,  vc  =  0,  hence  ns  —  0  for  all  J* s.  ns(J) 
is  plotted  in  Fig.  3  for  four  different  temperatures.  Here  J  is  normalized  with  the  Jc  at 
T  =  0.  At  T  =  0,  n3  =  ng  for  J  <  Jc(0).  This  is  what  has  been  predicted  by  the  two-fluid 
model. 

The  dependence  of  na  on  J  has  also  been  studied  in  the  Ginzburg-Landau  (GL)  theory 
[4,7].  Figure  4  shows  the  n3(J)/ns(J  =  0)  function  derived  from  the  GL  theory  [4]  with 
comparison  to  our  model.  In  our  model,  n3  decreases  smoothly  to  zero.  The  GL  theory 
predicts  that  there  is  a  critical  Jc  at  all  temperatures.  When  J  >  Jc(T),  ns  =  0,  or  n3 
becomes  an  imaginary  value.  This  n3(J)  function  was  derived  from  minimizing  the  Gibbs 
free  energy  g3.  When  minimizing  g3  with  respect  to  the  order  parameter  tf>  (n3  =  tftf*), 
only  a  first  few  terms  were  taken  in  expanding  g3  in  terms  of  power  series  of  rj>. 
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The  GL  theory  is  similar  to  the  macroscopic  quantum  model  (MQM)  [4]  for  homogeneous 
samples.  In  the  theory,  all  electrons  are  described  by  a  single  macroscopic  wavefunction 
<5.  Therefore  all  electrons  have  the  same  velocity,  same  as  v^.  The  distribution  n(v)  is 
taken  as  uqS(v  —  v^).  Hence  it  predicts  an  abrupt  drop  of  ns  to  zero  when  J  >  Jc(T)  for 
all  temperature  T’s.  This  is  not  what  has  been  observed  in  experiments.  Smooth  variation 
of  ns(J)  is  found  for  all  J  values.  If  the  wavefunction  of  each  electron  V’i  is  taken  into 
account  for  a  distribution  n(ut-),  where  =  <  ipi\v\tl>{  >,  a  modified  formula  can  be  derived 
from  statistical  physics. 

The  original  GL  theory  is  for  T  close  to  Tc ■  The  temperature  dependence  of  n3  assumed 
in  the  two-fluid  model  is  that 

n.(T)  =  no[l  -  (j~)4].  (15) 

The  ns(T)  function  is  compared  with  our  distribution  model  in  Fig.  5.  From  eq.  (14),  for 
a  fixed  J,  the  ns(T)  curve  matches  the  two-fluid  model  for  a  =  3/2. 


V.  Conductivities  <rs(J)  and  <rn(J) 


To  derive  the  constitutive  relations  of  superconductors  for  electromagnetic  fields,  we  will 
use  the  London’s  equations.  The  London’s  equations  are  derived  from  the  fundamental 
Newtonian  dynamics  and  the  Meissner  effect  [4].  They  can  also  be  derived  from  quantum 
mechanics  by  introducing  a  canonical  momentum  [3,8].  If  we  do  not  consider  the  Lorentz 
force  due  to  the  magnetic  field,  the  linear  London’s  equations  are  valid  [4]  from  the  Drude 
model  and  the  Newton’s  second  law. 


dJs  =  E 
8t  A 

„  T  B 
VxJs  =  -x 


(16) 

(17) 
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where  E  and  B  axe  the  total  electric  and  magnetic  fields  respectively,  Ja  is  the  current 
density  due  to  superelectrons,  and  A  =  rns/q2ns-  The  subscript  s  denotes  that  the  quantity 
is  of  superelectrons.  The  nonlinearity  will  be  included  in  A(J)  =  /xqA^,  and 


A(J)  = 


ms 


(18) 


q$ns(J)no 

Once  ns(J)  is  known,  A (J)  can  be  derived.  Hence,  the  nonlinear  constitutive  relations  are 
obtained.  Nonlinear  effects  come  in  A  and  ns. 

Substituting  eq.  (14)  of  ns(J)  in  A,  we  derive  the  A (J)  function,  which  is  plotted  in  Fig. 
6.  The  penetration  depth  goes  to  infinity  when  J  »  JC(T).  Conductivities  o-j(J)  and 
<rn(J)  are  derived  from  the  London  equations  for  time-harmonic  (e-w,;*)  fields 


and 


where 


<rn  = 


iq$ns 

(19) 

ums 

q2nnr 
m(l  —  twr) 

(20) 

-  />*•  -  C-ys-*--"™"”"*' 

and  the  number  density  of  normal  electrons  nn(J)  =  nQ  —  ns(J)  from  conservation  of 
charge.  The  above  equations  are  the  key  results  of  this  paper.  The  conductivities  axe 
plotted  in  Fig.  7a  and  7b,  where  r  is  the  transport  time  or  the  mean  scattering  time  of 
normal  electrons.  <ra(  J)  and  <rn(J )  reflect  the  current  dependence  of  na  and  nn  respectively. 
These  plots  axe  for  T  =  88K  very  close  to  the  critical  temperature  Tc  =  90  K.  The  nonlinear 
behavior  is  easy  to  see.  Total  conductivity  is  a  complex  number 

dr(J)  =  <rn(J)  +  <ra{  J)  (21) 


VI.  Geometry  Effect  in  Nonlinear  Relations 
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1.  Thin  wire 


For  a  thin  wire  superconductor  with  radius  a  «  A,  the  current  density  J  can  be  assumed 

independent  of  the  radius  p.  The  thin  wire  is  an  ideal  case  since  a  <<  A  is  not  practical. 

Usually  r  «  10-13  sec,  and  u>  <  10®  Hz,  hence  ut  «  1  which  is  the  quasi-static  case. 
2 

<7-n  ~  *  ^T.  Here  we  can  not  assume  |<Tj|  >>  |<rn|  since  we  are  studying  the  transition 
where  n3  may  become  very  small. 


Resistance  per  unit  length  of  the  wire  can  be  written  as 

R  =  [Kp  +  M  w  (22) 

The  kinetic  inductance  L ^  per  unit  length  of  the  wire  can  be  written  as 

Lk  =  +  (23) 

In  this  case,  the  internal  and  external  inductances,  £t-n  and  Lex,  related  to  the  energy 

stored  in  the  magnetic  'tel.*  axe  independent  of  7,  therefore  are  linear.  Figure  8  show  the 

calculated  R  —  7  and  Lj.  —  I  curves  from  eqs.  (22)  and  (23),  where  ur  =  10-3, 7  =  Jira?. 

The  R  —  J  curve  appears  very  nonlinear  because  T  =  88K  is  very  close  to  Tc  =  90K.  One 

interesting  behavior  in  the  L^  —  J  curve  is  that  a  maximum  Ls  appears  near  Jc.  This  can 

be  explained  from  eq.  (23).  At  J  «  Jc{T)t  ns  »  nnu>T,  hence  Ly.  ~  l/n5  is  small.  At 

J  >>  JC(T),  rinwr  »  n3,  hence,  Ly.  «  n5/(nn)3  is  also  small.  When  nnu>T  and  ns  are 

comparable,  a  maximum  of  Ly.  may  be  achieved. 


A  special  case  is  that  at  DC  where  u>  =  0.  From  eq.  (19),  crs  — »  oo.  From  the  first  London 
equation,  E  =  — =  0.  When  J  »  Jet  n3  — *  0,  then  it  is  possible  for  E  ^  0.  If  Js  has 
a  small  time  fluctuation,  then  ^  is  not  exactly  zero,  the  above  quasi-static  discussion  can 
be  applied. 


2.  Slab 
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The  above  discussion  has  not  taken  into  account  the  geometry  of  the  superconductor. 
For  finite  dimension  superconductors,  for  example,  a  slab  in  Fig.  9,  the  non-uniform 
distribution  of  current  density  Jy(x)  will  cause  the  nonlinearity  occur  earlier  in  the  R  —  I 
and  Lfc—I  curves.  The  interned  and  external  inductances  L^n  and  Lex  will  also  be  nonlinear 
since  Jy(x)  is  determined  by  /, 


J  dz  J  dxJy(x)  =  J  dz  J  dxd(x)E(x) 


(24) 


where  E(x)  is  the  electric  field.  For  a  uniform  current  density,  Jy(x)  =  JA,  J A  =  I /Azd 
where  d  is  the  thickness  of  the  slab.  For  non-uniform  Jy(z),  at  the  edge  of  the  slab,  Jy 
is  bigger  than  JA,  n3(J)  is  smaller  at  the  edge.  Therefore,  for  the  same  magnitude  of  I, 
inhomogeneous  Jy(x)  will  exceed  Jc  at  some  z’s  even  when  JA  is  less  than  Jc. 


If  a  three-dimensional  velocity  distribution  is  considered  for  a  current  flow  in  the  y  direc¬ 
tion, 

»(*.,«»,*«)  -  n0(5^?)Ie-”,K2+(«»-»Jl)2+»i]/2tBr  (25) 

Jy  =  qnQVjL  and 


limT-+Qn(v)  =  nQS(vx)S(vy  -  vA)S(vz) 


(26) 


fVc  tJvl~vy  t\!vc~vy~vx 

ns(Jy)=  dvy  /—-<&*/  /  -  n(v)dvz  (27) 

J-Vc 


with  vA  =  Jy/qnQ. 


Under  the  quasi-static  condition  u»r  <<  1,  if  a  current  I  is  applied  along  y  direction  in 
the  slab,  a  magnetic  field  Ez(x)  will  result  from  the  applied  current  Jy(x).  Note  that  in 
the  finite-width  slab  of  a  linear  superconductor,  even  at  w  =  0,  current  Jy(x)  has  a  non- 
uniform  distribution,  due  to  the  penetration  depth  A.  From  the  second  London  equation 
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(28) 


(17)  and  VxA  =  BtVxH  =  J,  under  the  London  gauge,  V  •  A  =  0  and 


A  =  -A  Js, 


V2(J,  A2(  J)]  =  Js  +  Jn  =  J  (29) 

where  A  is  the  vector  potential.  In  general,  this  equation  can  not  be  solved  in  a  closed-form. 


Let’s  first  consider  the  solution  of  a  linear  superconductor  for  ns  »  nnwr. 

_  .  .  I  cosh(z/A)  ,  .  . 

J»(X)  =  2Adnh(<l/2A)  }ar  WS,£/2 


(31) 

An  iterative  scheme  is  to  use  the  linear  solution  Jy(x)  in  eq.  (30)  to  calculate  ns(J)  in 
eq.  (14)  and  A(x)  in  eq.  (18).  Then  A  in  eq.  (30)  is  substituted  by  the  obtained  A(jc) 
to  calculate  Jy{x).  Then  eq.  (14)  is  used  to  calculate  n5[Jy(z)]  again.  This  procedure  is 
repeated  until  Jy(x)  converges.  Figure  10  shows  the  results  of  J(x)  and  <Ta(x)  obtained 
from  the  iterative  method. 


The  impedance  of  the  slab  is  defined  as 

%  = 

'd2/2dxa(x) 

and  resistance  R  =  Re(Z)  and  L j.  =  —Im(Z)/u>.  Figure  11  shows  the  R  —  J,  Ly.—  I 
relations  from  final  convergent  Jy(x). 

For  different  thicknesses  d,  the  nonlinear  curves  are  different.  We  have  found  that  for 
smaller  dj A,  R  —  I  and  —  I  relations  are  more  nonlinear.  This  is  understandable  since 
Jy  is  bigger  for  smaller  d  at  a  given  applied  current  I. 
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Surface  impedance  of  a  superconductor  is  defined  as 


=  Jf 


(33) 


where 


e  =  e  + 


id-(J) 


OJ 


(34) 


Surface  resistance  is  Ra  =  Re(Zs),  and  surface  inductance  is  La  =  —Im(Zs)/ iv.  For 
wr  =  10-3,  R$  and  La  axe  plotted  in  Fig.  12.  The  nonlinear  region  appears  near  Jc 
where  the  sample  is  partially  superconducting.  In  the  normal  or  superconducting  states, 
the  relations  are  linear.  The  La  is  different  at  normal  region  and  superconducting  region. 


Experimental  results  [13]  are  also  compared  in  Fig.  12c.  As  it  shows,  this  model  agrees 
with  the  trend  of  the  experimental  results.  Since  some  parameters  (a,  r,  etc.)  are  material 
dependent,  by  adjusting  these  parameters,  a  better  fit  between  the  theory  and  measure¬ 
ments  may  be  found.  Surface  impedance  Za  is  a  good  description  of  superconductors  since 
most  of  the  current  and  field  are  confined  within  the  penetration  depth  A  from  the  surface. 


Another  case  is  at  very  high  frequencies  and  crn  for  T  >  Tc  is  not  very  big  (which  is  true 

for  the  ceramic  high  Tc  superconductors)  so  that  e  can  not  be  neglected  in  eq.  (34)  for  Jy 
ft 

near  Jc .  Since  ^  0  in  the  Maxwell’s  equations,  a  wave  equation  has  to  be  considered. 
We  will  solve  the  guided  wave  case  where  the  electric  field  Ey  is  decaying  away  outside  the 
slab. 


The  wave  equation  for  E  is 

V2E  +  k$erE  =  0 


where 


fir 


=  1  + 


(35) 


(36) 
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If  we  replace  E  by  aJy  in  eq.  (35),  eq.  (35)  will  look  very  similar  to  eq.  (29).  We  can  also 
use  an  iterative  scheme  to  solve  for  the  nonlinear  cr  case.  First  we  assume  a  is  independent 
of  J,  solve  for  Jy(x).  Second  we  calculate  a(x)  =  o-[Jy(x)].  Third  we  solve  the  wave 
equation  for  the  inhomogeneous  medium  problem,  get  Jy(x).  These  steps  will  be  repeated 
until  Jy(x)  and  <r(x)  converge.  This  is  a  feed-back  process,  at  the  edge,  Jy  increases  will 
cause  ns  decreases,  hence  <rs  decreases,  hence  Jy  decreases,  until  stabilized. 

This  procedure  is  similar  to  solving  the  coupled  two  differential  equations  in  the  GL  theory. 
Where  the  equation  for  A  is  essentially  the  same  as  the  wave  equation  (35)  for  E  or  Jy.  The 
difference  is  the  second  equation  for  ip,  where  ip(A)  (n3(J)  )  will  be  derived.  We  derived 
ns(J)  from  the  Maxwellian  distribution,  while  the  GL  theory  derived  from  minimizing 
Gibbs  free  energy  at  equilibrium.  A  macroscopic  wave  function  (order  parameter)  is  used 
to  describe  all  electrons  in  the  GL  theory,  hence  the  distribution  of  velocity  diversity  is 
neglected. 

3.  Thin  film  strip 

The  slab  geometry  results  can  be  used  to  study  a  microstrip  geometry.  If  the  dimension 
of  z  in  Fig.9  is  squeezed  to  d  <<  A,  the  current  will  still  be  uniform  in  the  z  direction. 
Therefore  the  same  results  can  be  applied. 

This  nonlinear  study  can  also  be  applied  to  the  magnetic  field  H  dependence  of  R  and  L. 
Once  the  relation  between  the  magnetic  field  H  and  the  current  J  is  determined,  Hc  and 
Jc  are  related.  The  above  discussion  can  be  directly  applied. 

Although  the  model  is  classical,  the  corresponding  quantum  statistical  distribution  can  be 
used  to  derive  the  velocity  (energy)  distribution.  Discrete  distribution  may  be  needed  if 
the  energy  is  quantized. 
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VII.  Conclusion 


A  macroscopic  model  is  established  for  nonlinear  constitutive  relations  in  superconductors. 
Maxwellian  distribution  of  electron  velocities  is  used  to  derive  the  dependence  of  superelec¬ 
tron  density  ns  on  applied  current  density  J.  Complex  a(J)  is  obtained.  The  geometry  of 
a  superconductor  will  introduce  non-uniform  distribution  of  the  current  density  J.  There¬ 
fore  nonlinearity  will  be  enhanced  at  the  surface  of  the  device.  By  using  this  macroscopic 
model,  a  general  scheme  of  solving  for  electromagnetic  properties  of  superconductors  has 
been  proposed. 
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Figure  Captions 


Fig.  1.  Maxwellian  distribution  of  electron  velocities. 

Fig.  2.  n(«)  functions  for  different  average  velocity  tu’s. 

Fig.  3.  n.(J)/no  at  different  temperatures.  In  this  example,  a  =  3/2,  and  Te  =  90K. 
Jc(0)  is  from  eq.  (13). 

Fig.  4.  Comparison  between  the  GL  theory  and  this  model  for  n,(J)/n,(J  =  0)  at 
T  =  88K  .  Here  a  =  3/2,  Te  =  90K. 

Fig.  5.  Comparison  between  the  two-fluid  model  and  this  model  for  n,(T) /no  at 
J  =  0.  Here  Te  =  90k. 

Fig.  6.  Penetration  depth  A  as  a  function  of  current  density  J  at  T  =  88K.  Here 
a  =  3/2,  Tc  =  90K. 

Fig.  7.  (a)  super-conductivity  a,  as  a  function  of  J  at  T  =  88K.  Here  a  =  3/2, 
Te  =  90K.  (b)  normal  conductivity  <rn  as  a  function  of  J  at  T  =  88K.  Here  a  =  3/2, 
Te  =  9oK. 

Fig.  8.  (a)  Resistance  of  a  thin  wire  R  as  a  function  of  J  at  T  =  88K.  Here  a  =  3/2, 
Te  =  90K.  (b)  Kinetic  inductance  Lj,  of  a  thin  wire  as  a  function  of  /  at  T  =  88K. 
Here  a  =  3/2,  Te  =  90K. 

Fig.  9.  A  superconducting  slab  with  a  thickness  A. 

Fig.  10.  (a)  Current  density  Jy(z)  distribution  in  a  slab,  (b)  super-conductivity  <r,(x) 
distribution  in  a  slab.  These  are  for  T  =  80K  and  a  =  3/2,  T«  =  90K,  d/A(J  =  0)  =  2. 

Fig.  11.  For  d/A(J  =  0)  =  O.S,  1,2,  (a)  resistance  of  a  slab,  A,  as  a  function  of  applied 
current  intensity  I  at  T  =  80K;  (b)  kinetic  inductance  Lu  of  a  slab  as  a  function  of 


applied  current  intensity  I  at  T  —  80K.  Here  a  =  3/2,  T«  =  90K. 


Fig.  12.  (a)  Surface  resistance  Rt  as  a  function  of  J  at  T  =  88K  and  (b)  surface 
inductance  I,  as  a  function  of  J  at  T  =  88K,  Here  a  =  3/2,  Tc  =  90K.  (c)  R,  as  a 
function  of  surface  current  density  K  at  T  =  77k.  The  circles  axe  the  measured  data 
from  reference  [14],  The  sample  is  a  YBCO  at  /  =  10.4GHz  with  Te  =  92K.  The  solid 
curve  is  the  prediction  from  this  model,  where  a  =  3/2,  and  wt  =  2.8X10“*.  The 
surface  current  density  is  calculated  from  J  with  a  thickness  of  A(J). 
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Fig.  9  A  superconducting  slab 
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ELECTROMAGNETIC  PROFILE  RECONSTRUCTION 
USING  THE  RICCATI  EQUATION  APPROACH 

Jiqing  Xia  *,  Arthur  K.  Jordan  +  ,  Jin  A.  Kong 
Department  of  Electrical  Engineering  and  Computer  Science 
Massachusetts  Institute  of  Technology,  Cambridge,  MA  02139 
Permanent  address:  Center  for  Advanced  Space  Sensing 
Naval  Research  Laboratory,  Washington  DC  20375 


The  Riccati  differential  equation  for  reflection  coefficients  in  one-dimensional 
inhomogeneous  media  is  applied  to  the  electromagnetic  inverse  scattering  prob¬ 
lem.  Two  types  of  the  Riccati  equation  in  literature,  Schelkunoff’s  and  Red- 
heffer’s,  are  derived  and  distinguished.  Based  on  inverting  Redheffer’s  Riccati 
equation,  both  linear  and  non-linear  inversion  formulae  are  proposed.  These 
renormalized  perturbation  formulae  reconstruct  the  dielectric  profile  from  the  re¬ 
flection  coefficient  at  the  surface  of  the  medium.  Existing  inversion  formulae  (in¬ 
cluding  the  Born  approximation)  which  were  obtained  from  the  Green’s  function 
approach  and  the  Gel’fand-Levitan-Marchenko  (GLM)  theory  are  now  derived 
from  the  Riccati  equation.  Four  inversion  schemes  based  on  inverting  the  Riccati 
equation,  linearized  Redheffer’s  Riccati  equation,  linearized  Schelkunoff’s  Riccati 
equation,  non-linear  Redheffer’s  Riccati  equation  and  non-linear  Schelkunoff’s 
Riccati  equation  approaches,  axe  used  to  invert  several  dielectric  profiles.  Com¬ 
parison  and  summary  of  these  methods  axe  given.  All  these  methods  give  higher 
order  inversion  results  than  the  first-order  Born  approximation.  The  inversion  is 
performed  on  band-limited  reflection  coefficients  in  frequency  domain.  An  inverse 
Liouville  transform  is  introduced  to  rigorously  recover  the  geometric  lengths  from 
the  stretched  coordinates  in  the  inversion  procedure. 


Electromagnetic  Inverse  Scattering 
In  Remote  Sensing 
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Inverting  the  physical  parameters  which  usually  appear  as  inhomogeneous  pro¬ 
files  from  measurements  has  been  a  challenging  task  for  the  remote  sensing  com¬ 
munity.  Due  to  complexity  of  the  media,  the  forward  problem  in  remote  sensing 
of  forest,  vegetation  and  sea  ice,  typically  in  random  media,  has  attracted  re¬ 
searchers’  efforts  in  the  past  two  decades.  Techniques  which  serve  the  ultimate 
goal  of  solving  the  inverse  problem  have  also  evolved  and  entered  into  a  practical 
stage. 

In  general,  the  inverse  methods  can  be  summarized  in  three  categories.  The 
first  one  is  a  data  base  approach.  Whenever  a  measurement  result  comes  in, 
it  will  be  compared  with  all  the  pre-stored  patterns  in  a  library  to  identify  the 
possible  targets.  In  industry,  it  is  called  data  interpretation  which  strongly  relies 
on  previous  experience.  The  accumulation  of  data  in  the  library  is  a  learning 
process.  The  advantage  of  this  method  is  that  the  unknowns  can  be  identified 
exactly  if  such  target  has  been  seen  before  so  that  the  cause  and  the  result  are 
directly  connected.  The  disadvantage  is  that  if  the  data  base  is  not  big  enough, 
the  target  can  not  be  identified,  or  if  the  data  base  is  too  big,  it  will  take  too  much 
time  to  search  for  and  match  the  right  pattern. 
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The  second  method  is  the  iterative  approach.  A  forward  model  of  deducing  re¬ 
sults  from  causes  has  been  established  before  the  inversion.  When  a  measurement 
is  obtained,  a  set  of  guessed  values  for  the  unknowns  will  be  put  in  the  froward 
model  as  the  possible  causes.  If  the  predicted  result  of  the  forward  model  does 
not  match  the  measurement,  an  adjusted  set  of  parameters  will  be  used  as  a  new 
guess.  The  whole  procedure  will  be  repeated  until  the  match  is  found.  During 
iterations,  how  to  correct  the  error  to  reach  fastest  the  convergent  value  is  very 
important.  This  procedure  is  also  called  optimization  or  error  minimization,  which 
by  itself  is  an  active  research  area.  Iterative  methods  are  particularly  useful  when 
the  solution  of  an  integral  equation  cam  not  be  found  in  an  explicit  form.  As  an  ex¬ 
ample  of  the  iterative  methods,  the  Born  approximation  approximates  the  forward 
model  in  a  linearized  form,  then  iterates  towards  convergence.  The  advantage  of 
the  iterative  approach  is  that  an  accurate  final  solution  can  be  found  at  a  relative 
fast  (compare  to  data  base  search)  speed.  The  method  is  also  flexible  for  different 
kinds  of  targets,  unlike  the  data  base  search  where  the  target  has  to  have  been 
seen  before.  The  disadvantage  is  that  the  convergence  of  iteration  is  not  always 
guaranteed.  There  are  usually  more  than  one  minimum  of  the  cost  function.  The 
initial  guess  is  sometimes  so  critical  to  assure  the  right  convergent  result  that  a 
priori  information  needs  to  be  used.  Non-uniqueness  is  a  serious  problem  in  all 
inversion  methods.  Sufficient  number  of  measurements  are  helpful  of  reducing  the 
degree  of  non-uniqueness.  Formulation  of  the  forward  model  for  the  final  inversion 
equation  determines  how  efficient  the  inversion  scheme  is,  which  sometimes  also 
affect  the  degree  of  non- uniqueness. 

A  recently  developed  inversion  method  referred  to  as  the  renormalized  Source- 
Type  Integral  Equation  (STIE)  approach  [1]  solves  the  integral  equation  derived 
from  the  Green’s  function  [2]  without  the  linear  approximation.  The  STIE  ap¬ 
proach  formulates  an  exact  forward  model,  hence  is  not  restricted  to  low  contrast 
profiles  or  weak  scattering.  The  STIE  method  has  been  applied  to  profile  inver¬ 
sion  problems  of  the  soil  moisture  [3]  and  oil  formation  in  boreholes  [4],  where  the 
medium  properties  are  described  by  the  permittivity  and  conductivity  distribu¬ 
tions.  The  unknown  profiles  are  inverted  from  the  electromagnetic  measurements 
at  remote  observation  points. 

The  third  category  of  inverse  methods  is  the  explicit  solution  approach.  When 
one  solves  the  inversion  equation,  closed-form  solutions  may  be  obtained  if  the 
problem  is  formulated  in  certain  ways  and  the  data  is  in  certain  forms.  Examples  of 
this  method  include  the  Gel’fand-Levitan-Marchenko  (GLM)  theory  when  rational 
reflection  coefficient  is  used  and  methods  of  inverting  the  Riccati  equation  [5]. 
The  advantage  is  obviously  that  the  convergence  is  guaranteed  and  the  speed  is 
excellent.  Unfortunately,  only  a  very  few  practical  problems  can  be  solved  by 
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using  this  method.  It  is  usually  used  as  a  check  for  other  inversion  methods  or  for 
mathematical  studies. 

In  this  paper,  we  will  review  the  above  inverse  methods  with  our  focus  on 
iterative  methods,  especially  the  STIE  approach.  The  STIE  approach  is  extended 
to  the  case  of  a  general  background  medium.  The  formulation  of  the  inversion 
equation  is  exact.  Numerical  optimization  methods  are  employed  in  the  solution 
of  the  inversion  equation.  The  inversion  equation  also  has  an  explicit  dependence 
on  the  unknowns  to  be  inverted  for;  this  allows  one  to  compute  the  derivatives 
(including  higher  orders)  of  the  response  with  respect  to  these  unknowns  in  a  closed 
form.  Therefore,  higher-order  convergence  can  be  achieved  without  increasing 
much  of  the  computational  time.  By  pre-storing  the  elements  of  the  inversion 
equation  (which  depend  only  on  the  background  medium  and  are  independent  of 
the  unknown  profile),  the  method  does  not  require  the  solution  to  the  full  forward 
problem  repeatedly  as  in  the  case  of  the  Distorted  Born  approach  and  is  therefore 
faster  in  implementation.  An  algorithm  is  formulated  for  simultaneously  inverting 
the  unknown  permittivity  and  conductivity  profiles.  As  an  example,  a  soil  moisture 
profile  is  inverted  from  measurements  above  the  ground. 
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